Supershell structure in trapped dilute Fermi gases 
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We show that a dilute harmonically trapped two-component gas of fermionic atoms with a weak 
repulsive interaction has a pronounced super-shell structure: the shell fillings due to the spherical 
harmonic trapping potential are modulated by a beat mode. This changes the "magic numbers" 
occurring between the beat nodes by half a period. The length and amplitude of this beating mode 
depend on the strength of the interaction. We give a simple interpretation of the beat structure in 
terms of a semiclassical trace formula for the symmetry breaking U(3) — > SO (3). 

PACS numbers: 03.75.Ss 05.30.Fk 
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In finite systems of fermions, quantum effects lead to 
bunching of energy levels resulting in shell structure. 
Well-known examples are the shell structures in atoms 
or nuclei 0, determining their chemical properties and 
stability. In fact, shell effects and the so-called "magic 
numbers" , corresponding to spherical shell closings, have 
been discovered in a variety of other finite fermion sys- 
tems. Metal clusters in which the delocalized valence 
electrons are bound in the field of the metallic ions 
or quantum dots in semiconductor heterostructures 
are famous examples. More recent experimental progress 
makes it possible to study yet another species of finite 
quantal systems: atomic gases, often weakly interacting, 
confined e.g. by an optical dipole trap 0. 

In this paper, we show that a harmonically trapped 
gas of fermionic atoms interacting by a weak repulsive 
two-body force may exhibit super-shell structure: the 
shell oscillations of the spherical harmonic oscillator are 
modulated by a beat structure, whereby the positions of 
the magic numbers are shifted by half a period between 
successive beats. We can explain this surprising result 
semiclassically by the interference of diameter and cir- 
cle orbits surviving the breaking of the U(3) symmetry 
of the harmonic oscillator by the leading anharmonicity 
term in the mean field. 

Similar super-shell structure was predicted for metallic 
clusters |5j , inspired by a semiclassical analysis of Balian 
and Bloch in terms of the periodic orbits in a spherical 
cavity and observed experimentally 0. Analogous 
ideas could be applied to the description of shell structure 
in transport properties of quantum wires 

Let us consider a dilute gas of fermionic atoms, con- 
fined by a spherical harmonic potential modeling an ex- 
ternal trap JJ| , interacting through a repulsive zero-range 
two-body potential. The many-body Hamiltonian is 
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where a is the s-wave scattering length. The value and 
the sign of a can be varied either by changing the type 



of atoms or by applying a magnetic field: Feshbach reso- 
nances |{j allow to tune the scattering length from a large 
positive to a large negative value. (Here, we focus on the 
repulsive case.) 

Due to the Pauli principle the 5 interaction only ap- 
plies to fermions of pairwise opposite spin. We consider 
a fully unpolarized two-component system with two spin 
states, so that the total particle density is composed 
of two different densities of equal magnitude, n(r) = 
n^(r) + 71-1- (r) = 2n'(r). In the weak-interaction regime, 
the interaction energy density is given by gnJ (r)n^ (r) = 
gn 2 (r)/A, where the coupling strength parameter g is in- 
troduced by 

g = Anh 2 a/m. (1) 
This leads to the single-particle Hartree-Fock equation 
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V^(r) = e^(r) , (2) 



where Vho is the harmonic oscillator (HO) trap potential. 

The dilutcness condition necessary to treat the interac- 
tion as a two-body process is that the interparticle spac- 
ing n -1 / 3 is much larger than the range of the interaction 



and that 



<C 1. This dimcnsionless parameter also 



limits the life time of the two-component atomic fermion 
gas due to dimer formation, which is a three-body pro- 
cess ^lj- To guarantee that this condition is fulfilled, we 
calculate the central density n(0) in the Thomas-Fermi 
approximation and plot level curves of log(a 3 n(0)) in a 
g-N landscape seen in Fig. ^ The diluteness condition 
is seen to be fulfilled for all considered combinations of 
particle numbers and interaction strengths. This also 
outrules the possibility of phase separation discussed in 
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Assuming spherical symmetry, Eq.J2J) reduces to its 
radial part 
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with U(f) — gn^ (r) + ^mLd 2 r 2 being the effective mean- 
field potential. Each state has a (21 + l)-fold angular 
momentum degeneracy. We solve Eq. self-consistently 
on a grid. The interaction term is updated (with some 
weight factors) in each iteration according to gn'(r) = 
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FIG. 1: (Color online)Level curves of iolog (a 3 n (0)) calcu- 
lated within the Thomas-Fermi approximation, g — Antfa/m 
(h=ui=m=l) . 

After convergence is obtained, the ground-state energy 
of the A-particle system is given by (Ep = Fermi energy) 

E tot (g,N)= ]T ei-g f n 2 (r)d 3 r. (4) 

In general, the ground-state energy as a function of N can 
be written as the sum of a smooth average part and an 
oscillating part, E tot = E av + E osc . The oscillating part, 
referred to as the shell-correction energy, or shell energy 
in short, reflects the quantized level spectrum {e^}. For a 
non-interacting Fermi gas in a spherically symmetric 3D 
harmonic trap, the leading-order term for the average 
energy is found in the Thomas-Fermi approximation to 
be [l2j E\l = (3N) 4 / 3 Huj/4:. For the repulsive interacting 
case, we find E av (g > 0) oc N a with a larger exponent 
a > 4/3. However, Eq. J2J with an interaction term 
linear in the density is only valid for moderate g values 
and in practice we are close to a = 4/3 (e.g., a w 1.35 for 
g = 2). Contrary to the non-interacting case, and also to 
self-saturating fermion systems (such as nuclei and metal 
clusters) with a nearly constant particle density, it is not 
possible here to obtain the smooth part of the energy by 
a simple expansion in volume, surface and higher-order 
terms. We therefore perform a numerical averaging of 
the energy over the particle number N in order to 
extract its oscillating part. 

In the non-interacting case (5=0) the shell energy E osc 
oscillates with a frequency 27T3 1 / 3 rj 9.06 as a function 
of N 1 / 3 and has a smoothly growing amplitude oc N 2 / 3 . 
This follows from the exact trace formula 01 f° r Eosc of 
the 3D harmonic oscillator, whose leading-order term is 



given by 

E% c ~ (3A0 § |J £ ^ cos (2vrfc (3JV)*) . (5) 
fe=i 

Hereby k is the repetition number of the primitive clas- 
sical periodic orbit of the system with action Sq(E) = 
2-kE/uj. The argument of the cosine function in Eq. 
is simply k times Sq(E)/H, taken at the Thomas-Fermi 
value of the Fermi energy E F (N) = (3Ny/ 3 huj. The 
gross-shell structure is governed by the lowest harmonic 
with k — 1. 

Switching on the interaction, this scenario changes. A 
beating modulation of the rapid oscillations is found. In 
Fig. |21 we show the shell energy versus N 1 / 3 for three 
values of the interaction strength, g=0.2, 0.4 and 2. A 
beating modulation of the amplitude of the shell energy, 
i.e., a super-shell structure, is clearly seen to appear for 
all cases. At small particle numbers and particularly 
for small g values, the shell energy is that of the non- 
interacting system, given by Eq. (J5J. For larger inter- 
action strengths the super-shell structure is more clearly 
seen, and several beating nodes appear for g=2. With 
increasing interaction strength the amplitude of the shell 
energy oscillations becomes smaller. For example, for 
particle numbers around 80 3 « 500000, the amplitude of 
the shell energy is about 40 hoj, which is only about 10 -6 
of the total ground-state energy. 

Through Fourier analysis of the calculated shell en- 
ergy, two frequencies are seen to smoothly appear with 
increasing g value around the HO frequency (9.06), see 
Fig. The exact values of the two frequencies depend 
on the range of particle numbers included in the analy- 
sis. The super-shell features appear when the contribu- 
tion to the effective potential from the interaction, gn\ 
is sufficiently large, i.e., at large values of g and N. We 
also observe that (almost) until the first super-node, i.e., 
TV 1 / 3 w 28 in Fig. the magic numbers agree with the 
HO ones (a = 0). Between the first two super- nodes, i.e., 
28 < A 1 ' 3 < 49 in Fig. g| the magic numbers for the 
interacting system are situated in the middle of two HO 
magic numbers, i.e., they appear at the maxima of the 
fast shell oscillations. Then, after the second super-node 
they roughly agree with the unperturbed HO ones again. 

In the following we outline a semiclassical interpreta- 
tion of these features The U(3) symmetry of the 
unperturbed HO system is broken by the term SU = gn^ 
in (J2J), resulting in the SO (3) symmetry of the interact- 
ing system. The shortest periodic orbits in this system 
are the pendulating diameter orbits and the circular or- 
bits with a radius corresponding to the minimum of the 
effective potential in @ including the centrifugal term. 
These two orbits lead to the observed supershell beating. 
The above symmetry breaking has so far not been dis- 
cussed in the semiclassical literature. In a perturbative 
approach 0] , it can be accounted for by a group average 
of the lowest-order action shift AS(o) brought about by 
the perturbation of the system: (e^ As ^)oeU(3)- Hereby 
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o is an element of the group U(3) characterizing a mem- 4-dimensional manifold CP 2 |l5j|. which for a perturba- 
ber of the unperturbed HO orbit family (ellipses or cir- tion 8U{r) — er 4 can be done analytically |13| . 
cles). For the average it is sufficient to integrate over the 
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FIG. 2: (Color online) The oscillating part of the ground state energy in units of hui as a function of N 1/B for g = 0.2 (blue), 0.4 
(red) and 2 (green). The two lower curves are displaced by 400?kj and 600/kj, respectively. The vertical dotted lines correspond 
to the HO magic numbers Nmag = M(M + 1)(M + 2)/3 for M = 1, 2,. . . 
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FIG. 3: (Color online) Fourier spectra of the shell energy for 
g — 0.2 (blue), 0.4 (red) and 2 (green). A peak splitting 
around uj = 9.06 is resolved in the middle and bottom panels. 
The second harmonics (k — 2) are seen around uj — 18. 
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In the perturbative regime (e <g; 1) we find the follow- 
ing perturbed trace formula: 



osc 



(AO 



2 A 00 

m uj v - ^ 



(-1) 



k r 



2£7T 3 ^ k 3 
k=l 



sin 



kS, 



— sm 



kSd 

h 



(6) 



where kS c i and kS c are the classical actions of the diam- 
eter and circle orbits, respectively. In the limit e —> 0, 
their difference goes as k(S c — Sd) — > kenEp(N)/m 2 u) 5 , 
so that © tends to the pure HO limit JSJl. Extracting 
e from a polynomial fit to the numerical potential U(r) 
in ®, one can qualitatively describe the beating of the 
shell energy E osc (N). With e of order ~5x 10~ 4 g, Eq. 
(|t)|) approximately reproduces the curves seen in Fig. [21 
up to the beginning of the second supershell. It explains, 
in particular, also the phase change in the position of the 
magic numbers N mag shown in Fig. 

To cover larger values of e (and N), we have developed 
|13| an analytical uniform trace formula for the potential 
U(r) = mui 2 r 2 /2 + er 4 , which contains the contributions 
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of the 2-fold degenerate families of diameter and circular 
orbits to all orders in e. This is analogous to uniform 
trace formulae obtained earlier for U(l) 0] and U(2) 
symmetry breaking [Tt| . 
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FIG. 4: (Color online) An enlarged part of Fig. 2 for g = 
2. The circles mark the harmonic oscillator magic numbers 

N> mag • 

The beat structure in E osc has some similarities with 
that found in nuclei and metal clusters 0- There are, 
however, two essential differences. 1. Those systems are 
self-saturating and have steep mean-field potentials that 
can be modeled by a spherical cavity [6( . The present sys- 
tem, in contrast, has a mean field with much smoother 
walls that are dominated at large distances by the con- 
fining harmonic potential. 2. The super-shells in the 
cavity model come from the interference of the shortest 
periodic orbit families with three-fold degeneracy, as is 
usual in spherical systems 0,0]. Here, however, the 
gross-shell structure comes from the diameter and cir- 
cle orbits which are only 2-fold degenerate, whereas the 
fully 3- fold degenerate families of tori with rational ratios 
lo t : lo v — n : m of radial and angular frequency only con- 
tribute to the finer quantum structures at higher energies. 
As described in [13j, they bifurcate from the circle orbit 
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with repetition numbers k > 3 and can be included in 
the trace formula using standard techniques 0, 0, |2(J . 

As mentioned above, in a Fermi gas of atoms with re- 
pulsive interaction (a > 0), atoms can be lost through 
three-body recombination events. Two atoms with op- 
posite spin form a molecule while the third takes up en- 
ergy. Having a low recombination rate, and thus a long 
life time of the system, is desired. Petrov made an es- 
timate of the loss rate of particles, n/n rzs 111 (na 3 ) 2 e/h 
where e is the average kinetic energy of atoms. Tak- 
ing e « a realistic energy scale in current experi- 
ments, we estimate the life time of atoms in the trap to 
be 10 _6 s, 10~ 3 s, l(T 2 s for g = 2, 0.4, 0.2 respectively, 
when the number of particles is so large that the first 
node of the super-shell is reached. Hence the life time is 
longer when g (or a) is smaller, reflecting that the loss 
rate is proportional to a 6 . The temperature regime of 
this super-shell structure is below OAfiK. 

In conclusion, we have seen that the shell structure of 
fermions with weak, repulsive interactions in a harmonic 
trap shows a pronounced beating pattern, with the single 
shell positions changing by half a period length between 
the different beat nodes. A Fourier analysis of the os- 
cillating shell-correction part of the Hartree-Fock energy 
shows clear peaks at two slightly different frequencies. 
This is interpreted semiclassically by the interference of 
the shortest periodic orbits generated by the breaking 
of the U(3) symmetry of the non-interacting HO sys- 
tem, which are the families of diameter and circle orbits, 
through a uniform trace formula given fully in [3 

and, 

in the perturbative limit, in Eq. |JBJ. 

For very weak interactions, the splittings of the highly 
degenerate HO levels have earlier been calculated pertur- 
batively within the WKB approximation |2l|. However, 
the perturbative results do not apply for the interaction 
strengths where super-shell structure appears visible. 
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